// ///////////////////////////// MIT License //////////////////////////////////// //
//                                                                                //
// Copyright (c) 2010 David Zsolt Manrique                                        //
//                    david.zsolt.manrique@gmail.com                              //
//                                                                                //
// Permission is hereby granted, free of charge, to any person obtaining a copy   //
// of this software and associated documentation files (the "Software"), to deal  //
// in the Software without restriction, including without limitation the rights   //
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      //
// copies of the Software, and to permit persons to whom the Software is          //
// furnished to do so, subject to the following conditions:                       //
//                                                                                //
// The above copyright notice and this permission notice shall be included in     //
// all copies or substantial portions of the Software.                            //
//                                                                                //
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     //
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       //
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    //
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         //
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  //
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN      //
// THE SOFTWARE.                                                                  //
//                                                                                //
// ////////////////////////////////////////////////////////////////////////////// //

#include <iostream>
#include <ostream>
#include <sstream>
#include <vector>
#include <algorithm>
#include <complex>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <fstream>
#include <string>
#include <algorithm>
#include <iomanip>
#include <vector>
#include <list>
#include <set>
#include <cstdarg>
#include <stdexcept>

#include <Eigen/Eigen>
#include <cmdioutil.h>


struct SH
{
    Eigen::MatrixXcd S,H;
    double ef;
    void resize(int r,int c) { S.resize(r,c); H.resize(r,c); }
    void symetrize() { S = 0.5*(S+S.adjoint()); H = 0.5*(H+H.adjoint()); }
    double antisym_norm() const { return 0.5*( (S-S.adjoint()).norm()+(H-H.adjoint()).norm() );  } 
};

bool approxequal(const SH & sh1,const SH & sh2, double tol = 1.0e-4)
{
    return ( ( (sh1.S-sh2.S).norm() < tol) && ( (sh1.H-sh2.H).norm() < tol) );
}

struct blockSH : SH
{
    std::vector<std::pair<int,int> > orbital_range;
    int n;
    
    void half_range(int ii,int jj,int & i1,int & j1,int & i2,int & j2) const
    {
        ii = ii < 0 ? ii + n : ii - 1;
        jj = jj < 0 ? jj + n : jj - 1;

        if(ii>jj) std::swap(ii,jj);
        
        i1 = ii + 1;
        j1 = i1 + (jj-ii+1)/2 - 1;

        i2 = j1 + 1;
        j2 = i2 + (jj-ii+1)/2 - 1;
    }
    
    int range_length(int ii,int jj)
    {
        ii = ii < 0 ? ii + n : ii - 1;
        jj = jj < 0 ? jj + n : jj - 1;

        if(ii>jj) std::swap(ii,jj);
        
        return jj-ii+1;
    }
    void get_part(int ii,int jj,SH & sh) const
    {
        ii = ii < 0 ? ii + n : ii - 1;
        jj = jj < 0 ? jj + n : jj - 1;

        if(ii>jj) std::swap(ii,jj);
        
        int i1 = orbital_range[ii].first;
        int n1 = orbital_range[jj].first+orbital_range[jj].second-i1;
        
        sh.S = S.block(i1,i1,n1,n1);
        sh.H = H.block(i1,i1,n1,n1);
    }

    void get_part_2x2(int ii,int jj,SH & sh11,SH & sh12,SH & sh21, SH & sh22) const
    {
        ii = ii < 0 ? ii + n : ii - 1;
        jj = jj < 0 ? jj + n : jj - 1;

        if(ii>jj) std::swap(ii,jj);
        
        int i1 = orbital_range[ii].first;
        int n1 = orbital_range[jj].first+orbital_range[jj].second-i1;
        
        sh11.S = S.block(i1,i1,n1/2,n1/2);
        sh11.H = H.block(i1,i1,n1/2,n1/2);

        sh12.S = S.block(i1,i1+n1/2,n1/2,n1/2);
        sh12.H = H.block(i1,i1+n1/2,n1/2,n1/2);

        sh21.S = S.block(i1+n1/2,i1,n1/2,n1/2);
        sh21.H = H.block(i1+n1/2,i1,n1/2,n1/2);

        sh22.S = S.block(i1+n1/2,i1+n1/2,n1/2,n1/2);
        sh22.H = H.block(i1+n1/2,i1+n1/2,n1/2,n1/2);
        
    }

    void get_part(int ii1,int jj1,int ii2,int jj2,SH & sh) const
    {
        ii1 = ii1 < 0 ? ii1 + n : ii1 - 1;
        jj1 = jj1 < 0 ? jj1 + n : jj1 - 1;

        ii2 = ii2 < 0 ? ii2 + n : ii2 - 1;
        jj2 = jj2 < 0 ? jj2 + n : jj2 - 1;

        if(ii1>jj1) std::swap(ii1,jj1);
        if(ii2>jj2) std::swap(ii2,jj2);
        
        int i1 = orbital_range[ii1].first;
        int n1 = orbital_range[jj1].first+orbital_range[jj1].second-i1;

        int i2 = orbital_range[ii2].first;
        int n2 = orbital_range[jj2].first+orbital_range[jj2].second-i2;
        
        sh.S = S.block(i1,i2,n1,n2);
        sh.H = H.block(i1,i2,n1,n2);
    }
    
    void get_atom(int I,int J,SH & sh) const
    {
        I = I < 0 ? I + n : I - 1;
        J = J < 0 ? J + n : J - 1;

        int i,j,ni,nj;

        i = orbital_range[I].first;
        ni = orbital_range[I].second;
        j = orbital_range[J].first;
        nj = orbital_range[J].second;
        sh.S = S.block(i,j,ni,nj);
        sh.H = H.block(i,j,ni,nj);
        
        sh.ef = ef;
    }

    void print_pattern(std::ostream & os,double tol = 1e-6) const
    {
        int rowcount = 0;
        for(std::vector<std::pair<int,int> >::const_iterator it = orbital_range.begin(); it != orbital_range.end(); it++) 
        {
            rowcount++; if(rowcount==10) { os << std::endl; rowcount = 0; }
            int colcount = 0;
            for(std::vector<std::pair<int,int> >::const_iterator jt = orbital_range.begin(); jt != orbital_range.end(); jt++)
            {
                colcount++; if(colcount==10) { os << ' '; colcount = 0; }
                
                int i,j,ni,nj;
                i = it->first;
                ni = it->second;
                j = jt->first;
                nj = jt->second;
                if(H.block(i,j,ni,nj).norm() > tol)  os << 'X'; else os << '.';
            }

            os << "      ";
            for(std::vector<std::pair<int,int> >::const_iterator jt = orbital_range.begin(); jt != orbital_range.end(); jt++)
            {
                int i,j,ni,nj;
                i = it->first;
                ni = it->second;
                j = jt->first;
                nj = jt->second;
                if(S.block(i,j,ni,nj).norm() > tol)  os << 'X'; else os << '.';
            }
            
            
            os << std::endl;
        }
        
        os << std::endl << std::endl;
    }
};

void load(std::string fn,blockSH & sh)
{
    std::complex<double> ii(0,1.0);
    
    std::vector<std::string> lines;
    io::load(fn,"#",lines);

    sh.ef = 0.0;
    std::vector<std::string> efermi_block;
    io::parse_block("<fermi-energy>","</fermi-energy>",lines,efermi_block);
    for(std::vector<std::string>::iterator it = efermi_block.begin(); it != efermi_block.end(); it++)
    {
        std::istringstream iss(*it);
        if (!(iss >> sh.ef)) throw std::runtime_error("Not able to parse fermi energy block in the input file!");
        break;
    }

    std::vector<std::string> orbital_block;
    if(!io::parse_block("<index-orbital>","</index-orbital>",lines,orbital_block)) throw std::runtime_error("No or inproper index-orbital block in the input file!");
    if(!orbital_block.size()) throw std::runtime_error("No or inproper index-orbital block in the input file!");
    
    sh.orbital_range.clear();
    int index,orbital;
    int inc = 0;
    int orbital_index = 0;
    for(std::vector<std::string>::iterator it = orbital_block.begin(); it != orbital_block.end(); it++)
    if(util::match("% %",*it,index,orbital))
    {
            std::pair<int,int> p(orbital_index,orbital);
            orbital_index += orbital;
            sh.orbital_range.push_back(p);
            if(++inc != index) throw std::runtime_error("Index has to be in order from 1 in the index-orbital block!");
    }
    int n = orbital_index;
    sh.H.resize(n,n); sh.H.setZero(); sh.S = sh.H; sh.S.setIdentity();
    sh.n = sh.orbital_range.size();

    double re,im;
    int i,j;
    std::vector<std::string> overlap_block;
    if(!io::parse_block("<overlap>","</overlap>",lines,overlap_block)) throw std::runtime_error("There is no overlap matrix block or it is inproper in the input file!");
    for(std::vector<std::string>::iterator it = overlap_block.begin(); it != overlap_block.end(); it++)
    if(util::match("% % % % ",*it,i,j,re,im)) { if((i<n+1)&&(j<n+1)&&(i>0)&&(j>0)) sh.S(i-1,j-1)=re+ii*im; else throw std::runtime_error("Inconsistent overlap matrix index with index orbital block!");  }
    else if(util::match("% % % ",*it,i,j,re)) { if((i<n+1)&&(j<n+1)&&(i>0)&&(j>0)) sh.S(i-1,j-1)=re+ii*0.0; else throw std::runtime_error("Inconsistent overlap matrix index with index orbital block!"); }
    
    std::vector<std::string> hamiltonian_block;
    if(!io::parse_block("<hamiltonian>","</hamiltonian>",lines,hamiltonian_block)) throw std::runtime_error("There is no hamiltonian block or it is inproper in the input file!");
    for(std::vector<std::string>::iterator it = hamiltonian_block.begin(); it != hamiltonian_block.end(); it++)
    if(util::match("% % % % ",*it,i,j,re,im)) { if((i<n+1)&&(j<n+1)&&(i>0)&&(j>0)) sh.H(i-1,j-1)=re+ii*im; else throw std::runtime_error("Inconsistent hamiltonian index with index orbital block!"); }
    else if(util::match("% % % ",*it,i,j,re)) { if((i<n+1)&&(j<n+1)&&(i>0)&&(j>0)) sh.H(i-1,j-1)=re+ii*0.0; else throw std::runtime_error("Inconsistent hamiltonian index with index orbital block!"); }
    if(!hamiltonian_block.size()) throw std::runtime_error("There is no hamiltonian block or it is inproper in the input file!");

}

void save(std::string fn,blockSH & sh)
{
    std::complex<double> ii(0,1.0);
    std::ofstream os(fn.c_str());
    
    os << "<fermi-energy>" << std::endl;
    os << sh.ef << std::endl;
    os << "</fermi-energy>" << std::endl;
    os << std::endl;
    os << "<index-orbital>" << std::endl;
    for(int i=0; i<sh.orbital_range.size(); i++)    
    os << i+1 << " "  << sh.orbital_range[i].second << std::endl;
    os << "</index-orbital>" << std::endl;
    os << std::endl;
    os << "<overlap>" << std::endl;
    for(int i=0;i<sh.S.rows();i++)
    for(int j=0;j<sh.S.cols();j++)
    if(std::abs(sh.S(i,j))> 1e-7 )
    {
        os << "   " << i+1 << "    " << j+1 << "     "<< std::scientific << std::real(sh.S(i,j)) << "    " << std::scientific << std::imag(sh.S(i,j)) << std::endl;
    }
    os << "</overlap>" << std::endl;
    os << std::endl;
    os << "<hamiltonian>" << std::endl;
    for(int i=0;i<sh.H.rows();i++)
    for(int j=0;j<sh.H.cols();j++)
    if(std::abs(sh.H(i,j))> 1e-7 )
    {
        os << "   " << i+1 << "    " << j+1 << "     "<< std::scientific << std::real(sh.H(i,j)) << "    " << std::scientific << std::imag(sh.H(i,j)) << std::endl;
    }
    os << "</hamiltonian>" << std::endl;

    os.close();
    
}


void load_gollum(std::string fn,blockSH & sh)
{
    std::ifstream is(fn.c_str());
    std::complex<double> ii(0,1.0);
    std::string dump;
    if(!(is >> dump >> dump >> dump >> dump >> dump >> dump >> dump)) throw std::runtime_error("Wrong gollum input format: Syntax error near the beginning!");
    if(!(is >> dump >> dump >> dump >> dump >> dump >> dump >> sh.ef )) throw std::runtime_error("Wrong gollum input format: Syntax error near the beginning!");
    if(!(is >> dump >> dump >> dump >> dump >> dump >> dump )) throw std::runtime_error("Wrong gollum input format: Syntax errour near the beginning!");
    int rows,cols; 
    if(!(is >> dump >> dump >> rows )) throw std::runtime_error("Wrong gollum input format: Syntax error near the beginning");
    if(!(is >> dump >> dump >> cols )) throw std::runtime_error("Wrong gollum input format: Syntax error near the beginning");
    int index;
    int pindex = 1;
    int pi = 0;
    int n = rows;    
    sh.orbital_range.clear();
    for(int j = 0; j < n; j++)
    {
        if(!(is >> index >> dump >> dump >> dump >> dump )) throw std::runtime_error("Wrong gollum input format: Syntax error in the orbital block!");
        if(index != pindex)
        {
            std::pair<int,int> p(pi,j-pi);
            sh.orbital_range.push_back(p);
            pindex = index;
            pi = j;
        }
    }
    std::pair<int,int> p(pi,n-pi);
    sh.orbital_range.push_back(p);
    sh.n = sh.orbital_range.size();

    if(!(is >> dump >> dump >> dump >> dump >> dump >> dump )) throw std::runtime_error("Wrong gollum input format: Syntax error near k-block!");
    //int rows,cols;
    if(!(is >> dump >> dump >> rows  )) throw std::runtime_error("Wrong gollum input format: Syntax error near k-block!");
    if(!(is >> dump >> dump >> cols  )) throw std::runtime_error("Wrong gollum input format: Syntax error near k-block!");
    for(int j = 0; j < rows; j++)
        if(!(is >> dump >> dump >> dump  )) throw std::runtime_error("Wrong gollum input format: Syntax error in k-block!");
    if(rows > 1) std::cerr << "Warning: Gollum input delivers multiple k point matrix. This load does not handle that, instead try separate xml inputs for each k-point." << std::endl;
    if(!(is >> dump >> dump >> dump >> dump >> dump >> dump  )) throw std::runtime_error("Wrong gollum input format: Syntax error near matrix block!");
    //int rows,cols;
    if(!(is >> dump >> dump >> rows  )) throw std::runtime_error("Wrong gollum input format: Syntax error near matrix block!");
    if(!(is >> dump >> dump >> cols  )) throw std::runtime_error("Wrong gollum input format: Syntax error near matrix block!");
    sh.H.resize(n,n); sh.H.setZero(); sh.S = sh.H;
    for(int k = 0; k < rows; k++)
    {
        int i,j;
        double res,ims,reh,imh;
        if(!(is >> dump >> i >> j >> res >> ims >> reh >> imh  )) throw std::runtime_error("Wrong gollum input format: Syntax error in matrix block!");
        sh.H(i-1,j-1)=reh+ii*imh;
        sh.S(i-1,j-1)=res+ii*ims;

    }
    is.close();
    
}


template<typename scattering>
void out_coef(double e, scattering & s,std::ostream & os,bool comment = false,bool advanced = false,int w = 18,int p = 10)
{
        //os.flags(std::ios::scientific);
        os.flags(std::ios::fixed);
        if(comment)
        {       
                std::ostringstream oss;  oss << "E-Ef(=" << std::setprecision(2)<< std::fixed << e <<")";
                os << '#';
                os << std::right << std::setw(w-1) << std::setprecision(p) << oss.str() << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << "Tran a->b" << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << "Tran b->a" << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << "Ref a->a" << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << "Ref b->b" << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << "a oc" << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << "b oc" << ' '; 
//                os << std::right << std::setw(w) << std::setprecision(p) << "Max tt' a->b" << ' '; 
//                os << std::right << std::setw(w) << std::setprecision(p) << "Max tt' b->a" << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << "max abs S*S'-I" << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << "eiva tt' a->b" << ' '; 
                os << std::endl;
        }
        else
        {   
            os << std::right << std::setw(w) << std::setprecision(p) << e << ' ';
            if(s.ano*s.bno !=0)
            {
                os << std::right << std::setw(w) << std::setprecision(p) << s.coefba << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << s.coefab << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << s.coefaa << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << s.coefbb << ' ';
            } 
            else
            {
                os << std::right << std::setw(w) << std::setprecision(p) << 0.0 << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << 0.0 << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << (double)s.ano << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << (double)s.bno << ' ';
                s.checknorm = 0.0;
            }
            os << std::right << std::setw(w) << std::setprecision(p) << s.ano << ' '; 
            os << std::right << std::setw(w) << std::setprecision(p) << s.bno << ' '; 

            if(s.ano*s.bno ==0) s.checknorm = 0.0;
            os << std::right << std::setw(w) << std::setprecision(p) << std::scientific << s.checknorm << ' '; 
            os << std::endl;                

/*          for(int i = s.ttabe.rows()-1; i >= 0; i--)
            os << std::right << std::setw(w) << std::setprecision(p) << std::scientific << s.ttabe(i).real() << ' ';  
            os << std::endl;                
            
            if(advanced)
            {
                os << std::endl;
                os << "advanced: scatterer wv for scatterer channles (experiment 2010jul16)" << std::endl;
                for(int i = s.ttabe.rows()-1; i >= 0; i--)
                os << std::right << std::setw(33) << std::setprecision(10) << std::scientific << s.ttabe(i).real() << ' ';  
                os << std::endl;                

                for(int i = 0; i < s.yml.rows(); i++)
                {
                    for(int j = s.yml.cols()-1; j >= 0; j--)
                    {
                        os << std::right << std::setw(15) << std::setprecision(5) << std::scientific << s.yml(i,j).real() << ' ';  
                        os << std::right << std::setw(17) << std::setprecision(5) << std::scientific << s.yml(i,j).imag() << ' ';  
                    }
                    os << std::endl;
                }
                os << std::endl;
                os << std::endl;
            }
*/

        }
        

}



template<typename lead,typename scatterer>
void get_amb(int argc, char* argv[],lead & a,scatterer & m,lead & b,double & emin_p, double & emax_p,double & ef_p,int & nume_p,double &tol_p,bool &advanced)
{
        cmd::arguments arg(argc,argv);
        cmd::switcharg<false> help(arg,"-h","--help", "print help message");
        cmd::switcharg<false> info(arg,"-i","--info", "print info/example message");

        cmd::vararg<double> emin(arg,"-emin","--emin","min of e-range (ef will be set to zero), <arg>=-1.0",-1.0);
        cmd::vararg<double> emax(arg,"-emax","--emax","max of e-range (ef will be set to zero), <arg>=1.0",1.0);
        cmd::vararg<int> nume(arg,"-n","--nume","number of energy intervals, <arg>=1",1);
        
        cmd::vararg<util::pair_expression<int,int,':'> > ha(arg,"-ha","--ha","lead-a uc pair atom index range, 1:3 or -1:-4","");
        cmd::vararg<util::pair_expression<int,int,':'> > hb(arg,"-hb","--hb","lead-b uc pair atom index range, 1:4 or -1:30","");

        cmd::vararg<util::pair_expression<int,int,':'> > ah(arg,"-ah","--ah","lead-a uc pair atom index range (h taken from opposite side)","");
        cmd::vararg<util::pair_expression<int,int,':'> > bh(arg,"-bh","--bh","lead-b uc pair atom index range (h taken from opposite side)","");

        cmd::varargxx<util::pair_expression<int,int,':'>,std::string> la(arg,"-la","--lead-a","replace i:j h range and lead-a ham. file","","");
        cmd::varargxx<util::pair_expression<int,int,':'>,std::string> lb(arg,"-lb","--lead-b","replace i:j h range lead-b ham. file","","");

        cmd::vararg<double> tol(arg,"-tol","--tol","channel separation tolerance def.: 1e-8",1.0e-8);

        cmd::switcharg<false> rand(arg,"-r","--rand", "default diagonal randomizer");        
        cmd::vararg<double> randa(arg,"-ra","--randa","randomize lead-a, to avoid singular coupling, <arg>=0.0 ",0.0);
        cmd::vararg<double> randb(arg,"-rb","--randb","randomize lead-b, to avoid singular coupling, <arg>=0.0",0.0);


        cmd::varargxx<util::pair_expression<int,int,':'>,double> homo(arg,"-homo","--homo","homo level shift of the selected range","",0.0);
        cmd::varargxx<util::pair_expression<int,int,':'>,double> lumo(arg,"-lumo","--lumo","lumo level shift of the selected range","",0.0);

        
        cmd::switcharg<false> gollum(arg,"-g","--gollum", "hamiltonian file in gollum format");
        cmd::switcharg<false> symmetrize(arg,"-s","--sym", "symetrize hamiltonian file before calculation");
        cmd::switcharg<false> adv(arg,"-a","--adv", "Advanced analysis and outputs");

        cmd::posarg<std::string> ham_fn(arg,"hamiltonain input file name", "");

        if(*help)
        {
                arg.print_help(std::cout,"transportab");
                std::cout << "examples:" << std::endl;
                std::cout << "transportab hamiltonian.dat -ha 3:4 -hb -4:-3 -n 100" << std::endl;
                std::cout << "transportab hamiltonian.dat -ha 3:4 -hb -4:-3 -n 100 -r" << std::endl;
                std::cout << "version 0.95b" << std::endl;
                exit(EXIT_SUCCESS);
        }
        if(*info)
        {
                std::cout << "version 0.95b" << std::endl;
                std::cout << "Copyright (c) 2010 David Zsolt Manrique (david.zsolt.manrique@gmail.com)"  << std::endl;
                
                //return EXIT_SUCCESS;
                exit(EXIT_SUCCESS);
        
        }                                                                                                                
        
        if(!~ham_fn) throw std::runtime_error("Give hamiltonian file!");
        
        blockSH sh;

        if(!*gollum) load(*ham_fn,sh); else load_gollum(*ham_fn,sh);
        
        if(*symmetrize) sh.symetrize();
        else if(sh.antisym_norm() > 1e-6) throw std::runtime_error("Hamiltonian or Overlap is not hermitian!"); /*check if gramm matrix e?*/

        blockSH sh_la,sh_lb;
        if(~la) { if(!*gollum) load(la.var2,sh_la); else load_gollum(la.var2,sh_la);  }
        if(~lb) { if(!*gollum) load(lb.var2,sh_lb); else load_gollum(lb.var2,sh_lb);  }

        //////////////////////////////////////////////////////////////////

        int ia,ja;
        if(~ha) { ia=std::get<0>(*ha); ja=std::get<1>(*ha); }
        else if(~ah) { ia=std::get<0>(*ah); ja=std::get<1>(*ah); }
        else if(~la) { ia=std::get<0>(la.var1); ja=std::get<1>(la.var1); }
        else throw std::runtime_error("Give -ha -ah range or -la range lead ham. file!");

        int ib,jb;
        if(~hb) { ib=std::get<0>(*hb); jb=std::get<1>(*hb); }
        else if(~bh) { ib=std::get<0>(*bh); jb=std::get<1>(*bh); }
        else if(~lb) { ib=std::get<0>(lb.var1); jb=std::get<1>(lb.var1); }
        else throw std::runtime_error("Give -hb -bh range or -lb range lead ham. file!");
        
        int im,jm;
        if(~la)
        {
            SH h11,h12,h21,h22; sh_la.get_part_2x2(ia,ja,h11,h12,h21,h22);

            int na = sh_la.range_length(ia,ja);
            
            SH h00;
            sh.get_part(1,na/2,h00);

            double de = (h00.H - h22.H).diagonal().array().sum().real() / h22.S.diagonal().array().sum().real();
            
            a.h = h22.H + de * h22.S;
            a.s = h22.S;
            a.v = h12.H + de * h12.S;
            a.t = h12.S;
            
            ia = 1;
            ja = na/2;
            im = ja+1;
        }
        else if(~ha)
        {
            SH h11,h12,h21,h22;
            sh.get_part_2x2(ia,ja,h11,h12,h21,h22);

            a.h = h11.H;
            a.s = h11.S;
            a.v = h12.H;
            a.t = h12.S;
            
            int ia2,ja2;
            sh.half_range(ia,ja,ia,ja,ia2,ja2);
            im = ja+1;
        }
        else if(~ah)
        {
            SH h11,h12,h21,h22;
            sh.get_part_2x2(ia,ja,h11,h12,h21,h22);

            a.h = h22.H;
            a.s = h22.S;
            a.v = h12.H;
            a.t = h12.S;

            int ia1,ja1;
            sh.half_range(ia,ja,ia1,ja1,ia,ja);
            im = ja+1;
        }


        if(~lb)
        {
            SH h11,h12,h21,h22; sh_lb.get_part_2x2(ib,jb,h11,h12,h21,h22);

            int nb = sh_lb.range_length(ib,jb);
            
            SH h00;
            sh.get_part(-1,-nb/2,h00);

            double de = (h00.H - h11.H).diagonal().array().sum().real() / h11.S.diagonal().array().sum().real();
            
            b.h = h11.H + de * h11.S;
            b.s = h11.S;
            b.v = h12.H + de * h12.S;
            b.t = h12.S;

            ib = -1;
            jb = -nb/2;
            jm = jb-1;
        }
        else if(~hb)
        {
            SH h11,h12,h21,h22;
            sh.get_part_2x2(ib,jb,h11,h12,h21,h22);

            b.h = h22.H;
            b.s = h22.S;
            b.v = h12.H;
            b.t = h12.S;

            int ib1,jb1;
            sh.half_range(ib,jb,ib1,jb1,ib,jb);
            jm = jb-1;
        }
        else if(~bh)
        {
            SH h11,h12,h21,h22;
            sh.get_part_2x2(ib,jb,h11,h12,h21,h22);

            b.h = h11.H;
            b.s = h11.S;
            b.v = h12.H;
            b.t = h12.S;

            int ib2,jb2;
            sh.half_range(ib,jb,ib,jb,ib2,jb2);
            jm = jb-1;
        }

        SH shm,shva,shvb;
        sh.get_part(im,jm,shm);
        sh.get_part(ia,ja,im,jm,shva);
        sh.get_part(im,jm,ib,jb,shvb);    
        
        m.h = shm.H;
        m.s = shm.S;
        
        m.va = shva.H;
        m.ta = shva.S;
        
        m.vb = shvb.H;
        m.tb = shvb.S;


        if(*rand)
        {
            Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> r = a.v;
            r.setRandom();
            a.v += 0.0001 * r.real().cast<std::complex<double> > ();
            r = b.v;
            r.setRandom();
            b.v += 0.0001 * r.real().cast<std::complex<double> > ();
        }
        
        if(~randa)
        {
            Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> r = a.v;
            r.setRandom();
            a.v += *randa * r.real().cast<std::complex<double> > ();
        }

        if(~randb)
        {
            Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> r = b.v;
            r.setRandom();
            b.v += *randb * r.real().cast<std::complex<double> > ();
        }

        tol_p = *tol;
        emin_p = *emin;
        emax_p = *emax;
        ef_p = sh.ef;
        nume_p = *nume;
        advanced = *adv;
}

/*
#include <iostream>
#include <ostream>
#include <sstream>
#include <vector>
#include <algorithm>
#include <complex>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <fstream>
#include <string>
#include <algorithm>
#include <iomanip>
#include <vector>
#include <list>
#include <set>
#include <cstdarg>
#include <stdexcept>

#include <Eigen/Eigen>

#include <cmdioutil.h>
*/

